rm(list=ls())
source("extras/packages.R")


## please download from https://www.worldvaluessurvey.org/WVSOnline.jsp
dat = readRDS("data/WVS_Cross-National_Wave_7_rds_v6_0.rds") 



### keep only relevant columns
# Q171: attend religious service
# Q289: religious denomination


cols = c("B_COUNTRY", "B_COUNTRY_ALPHA", "Q171", "Q289")
dat = dat[,cols]
colnames(dat) = c("code", "alpha", "relservice", "relgroup")


### recode vars to ensure correct values
dat$relgroup[!(dat$relgroup %in% 0:9)] = 0
dat$relgroup2pt = car::recode(dat$relgroup, "0=0 ; 1:9 = 1")

dat$relservice[!(dat$relservice %in% 1:7)] = 0
#dat$relservice[dat$relgroup2pt == 0] = NA
dat$relservice2pt = car::recode(dat$relservice, "1:3 = 1; 4:7 = 0")



### collapse data
temp <- dat %>%
  group_by(code, alpha) %>%
  dplyr::summarise(relservice2pt = round(mean(relservice2pt, na.rm=TRUE)*100,2),
                   relgroup2pt = round(mean(relgroup2pt, na.rm=TRUE)*100,2))


### create country codes
## please download from https://www.worldvaluessurvey.org/WVSDocumentationWVL.jsp
ref = read.xlsx2("data/F00012255-WVS_TimeSeries_1981_2020_CountrySpecificCodes.xlsx", sheetIndex = 1,
                 header = FALSE)
colnames(ref) = c("code", "region")


### merge country names and collapsed data
temp = merge(x = temp, y = ref,
             all.x = TRUE, all.y = FALSE)


### sort
temp = temp[order(temp$relgroup2pt, decreasing = TRUE),]

temp = temp[1:30, ]


### create long format
ctydat = data.frame(country = c(temp$region, temp$region), 
                    var = rep(c("Identifying with a Religion", 
                                "Attending Service at least Monthly"), 
                              each = nrow(temp)), 
                    value = c(temp$relgroup2pt, temp$relservice2pt))

ctydat$country = factor(ctydat$country, levels = temp$region[nrow(temp):1])

ctydat_relgroup = subset(ctydat, var == "Identifying with a Religion")
ctydat_relservice = subset(ctydat, var == "Attending Service at least Monthly")


p <- ggplot(ctydat) + 
  geom_segment(data = ctydat_relservice,
               aes(x = value, xend = ctydat_relgroup$value, 
                   y = country, yend = ctydat_relgroup$country), 
               color = "#aeb6bf",
               linewidth = 1.5,
               alpha = .5) + 
  geom_point(aes(x = value, y = country, color = var, shape = var), size = 2.5, show.legend = TRUE) +
  scale_x_continuous(name = "Percentage",
                     breaks = seq(0, 100, 10),
                     limits = c(0, 100)) + 
  labs(color = "",
       shape = "",
       y = "Country") +
  theme_bw() +
  theme(legend.position = "bottom")

ggsave(p, filename="Figure-1.pdf", height = 8, width = 9)

